# === Climate Network Replication Script ===
# Reproduces figures and tables for: Dynamic Networks of Negotiation for International Climate Change Cooperation
# Authors: Zack Almquist, Benjamin Bagozzi, Daria Blinova, Zach Brown

# === 1. Libraries ===
library(tidyverse)
library(network)
library(sna)
library(rnaturalearth)
library(rnaturalearthdata)
library(sp)
library(rgdal)
library(sf)
library(reshape2)
library(gridExtra)
library(cowplot)
library(ergm)
library(dnr)
library(ggplot2)
library(arm)


# === 2. Load and Prepare Data ===
setwd("Your_Directory_Here")
big_data <- read.csv("unfccc_symmetric_z.csv")
#Error for Western Samoa
big_data <- big_data %>% filter(country1 != "Western Samoa", country2 != "Western Samoa")
covar<-read.csv("UNFCCCcountry_covariates.csv")

# === 3. Build Annual Networks ===
edge_list <- lapply(unique(big_data$year), function(year) {
  x <- big_data[big_data$year == year, ]
  nodes <- unique(unlist(c(x[,70], x[,75])))
  m1 <- match(x[,70], nodes)
  m2 <- match(x[,75], nodes)
  mentions <- x$mentions
  df <- data.frame(from = as.numeric(m1), to = as.numeric(m2), weight = mentions)
  attr(df, "vertex.names") <- nodes
  df
})

net_cop <- lapply(edge_list, function(x) {
  net <- network.initialize(length(attr(x, "vertex.names")), directed = FALSE)
  network.vertex.names(net) <- attr(x, "vertex.names")
  apply(x[x$weight != 0, ], 1, function(row) net[row[1], row[2]] <<- 1)
  net
})

all_edges <- do.call(rbind, edge_list)
cumulative_edges <- aggregate(weight ~ from + to, data = all_edges, sum)
cumulative_net <- network.initialize(length(attr(edge_list[[1]], "vertex.names")), directed = FALSE)
network.vertex.names(cumulative_net) <- attr(edge_list[[1]], "vertex.names")
for (i in 1:nrow(cumulative_edges)) {
  cumulative_net[cumulative_edges$from[i], cumulative_edges$to[i]] <- cumulative_edges$weight[i]
}
#Note the two different bins here: net_cop contains the cumulative edges, net_cop_og contains just the annual networks
net_cop_og<-net_cop
net_cop[[length(net_cop) + 1]] <- cumulative_net

# === 4. Visualizations ===

# === Figure 1: World Map of Climate Network ===
world <- ne_countries('large')
nations_df<-big_data[,70:66]
nations_df<-unique(nations_df)

world_points <- data.frame(country=nations_df$country1,lat=nations_df$country1_lat,lon=nations_df$country1_lng)
world_points <- rbind(world_points, list('Samoa', 13.7590,172.1046))
sp<-SpatialPointsDataFrame(world_points[,3:2],
                           proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
                           data=world_points)

sp_sf <- st_as_sf(sp)
sp_sf <- st_transform(sp_sf, st_crs(world))

# Get coordinates from 'sp_sf'
coord <- st_coordinates(sp_sf)

# Set COP
i <- 14

# Open PDF for saving the plot
pdf("SpatialNetwork.pdf", width = 20, height = 20)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0) + 0.1)

# Plot the 'world' object (assumed to be 'sf')
plot(st_geometry(world), border = rgb(0, 0, 0, .4)) # Use st_geometry for sf plotting

# Plot the network on top of the world map
plot(net_cop[[i]], coord = coord, displayisolates = TRUE,
     edge.col = rgb(0, 0, 0, .06),
     edge.lwd = .8,
     vertex.cex = .1,
     vertex.col = rgb(1, 0, 0, .8),
     vertex.border = rgb(0, 0, 0, .3), jitter = TRUE,
     label = network.vertex.names(net_cop[[i]]),
     label.pos = 4,
     label.cex = .6,
     label.col = rgb(0, 0, 0, .9),
     new = FALSE
)

# Close the PDF device
dev.off()

# === Figure 3: Annual Network Snapshots ===
par(mfrow=c(1,1))
i<-1
coord<-plot(net_cop_og[[i]])
year<- 2010:2019
year<-append(year, 2021:2023)
cop<-16:28

pdf("Annual_Networks.pdf")
par(mfrow=c(3,5), mar = c(0, 0, 0, 0) + 0.1)
for(i in 1:(length(net_cop_og))){
  plot(net_cop_og[[i]],coord=coord,displayisolates=FALSE,edge.col=rgb(0,0,0,.1))
  title(sub = paste("COP ",cop[i]," (",year[i],")",sep=""),line=-1)
  box()
}
dev.off()
dev.off()

# === Figure 2: Matrix Plot ===
p<-vector("list",length(cop))
for(i in 1:length(p)){
  e<- sna::evcent(net_cop[[i]])
  o<-order(e,decreasing = TRUE)
  longData<-melt(net_cop[[i]][o,o])
  p[[i]]<-(ggplot(longData, aes(x = Var2, y = Var1)) + 
             geom_raster(aes(fill=value)) + 
             scale_fill_gradient(low="grey90", high="red") +
             labs(x="", y="", title="Matrix Plot") +
             theme_bw() + 
             theme(legend.position = "none")+
             theme(axis.text.x=element_text(size=3, angle=90, vjust=-.1),
                   axis.text.y=element_text(size=3),
                   plot.title=element_text(size=11))+
             ggtitle(paste("COP",cop[i]))
  )
  print(p[[i]])
}

cowplot::plot_grid(plotlist = p, ncol = 5)

ggsave(file = "Matrix_Plots.pdf", arrangeGrob(grobs = p, ncol = 5),dpi=600,width = 24, height = 16)

# === 5. Descriptive Statistics (Table 1) ===

results2 <- data.frame(Graph = 1:length(net_cop_og), Clustering_Coefficient = NA, SE_CC = NA, Density = NA, SE_Den = NA, Mean_Deg = NA, SE_MDeg = NA)

for (i in seq_along(net_cop_og)) {
  graph <- net_cop_og[[i]]
  
  # Clustering coefficient
  cc <- gtrans(graph, mode = "graph")
  
  # Number of nodes
  n <- network.size(graph)
  
  # SE for clustering coefficient
  se_cc <- sqrt((cc * (1-cc)) / n)
  
  # Density
  den<-gden(graph, mode = "graph")
  
  # SE for Den
  se_den<-sqrt((den * (1-den)) / n)
  
  #Mean Deg
  degrees <- degree(graph, gmode = "graph")
  n2 <- length(degrees)
  mean_deg<-mean(degrees)
  se_deg <- sd(degrees) / sqrt(n2)
  
  # Store results
  results2[i, "Clustering_Coefficient"] <- cc
  results2[i, "SE_CC"] <- se_cc
  results2[i, "Density"]<-den
  results2[i, "SE_Den"]<-se_den
  results2[i, "Mean_Deg"]<-mean_deg
  results2[i, "SE_MDeg"]<-se_deg
}

write.csv(results2, "Descriptives_Table.csv")

# === 6. DNR Modeling (Table 2) ===

covar_t <- big_data %>%
  select(
    country = country1,
    year,
    lagGDP = c1wdi_loggdppc_lag,
    lagpop = c1wdi_logpop_lag,
    #lagpopden = lagpopden1,
    lagGDPgrowth = c1wdi_gdpgr_lag
  )
covar_t <- covar_t %>%
  arrange(country, year)

covar_t <- covar_t %>%
  distinct(country, year, .keep_all = TRUE)

world_points <- data.frame(country=covar$country,lat=covar$lat,lon=covar$lng)
sp<-SpatialPointsDataFrame(world_points[,3:2],
                           proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
                           data=world_points)

## Distance matrix in km
d<-spDists(sp,longlat = TRUE)
rownames(d)<-colnames(d)<-covar$country

model.terms <- c("triadcensus.003", "triadcensus.012",
                 "triadcensus.102", "triadcensus.021D")
model.formula <- net~triadcensus(0:3)-1
graph_mode <- 'digraph'
group <- 'dnc'
alpha.glmnet <- 1
directed <- TRUE
method <- 'bayesglm'
maxlag <- 3
lambda <- NA
intercept <- NA
cdim <- length(model.terms)

#Need to run as a chunk from here to get paper values, includes the seed**************************************
set.seed(1982)
lagmat <- matrix(sample(c(0,1),(maxlag+1)*cdim, replace = TRUE),ncol = cdim)
ylag <- rep(1,maxlag)
exvar <- NA

tmp<-net_cop_og
tmp<-lapply(tmp,function(x){
  covar_t[is.na(covar_t)]=0
  x%v%"g77china"<-covar$g77china
  x%v%"eu"<-covar$eu
  x%n%"log_d"<-log(d+1)
  x%v%"african_group_negotiators"<-covar$african_group_negotiators
  x%v%"lagGDPpc"<-covar_t$lagGDP
  x%v%"lagpop"<-covar_t$lagpop
  #x%v%"lagpopden"<-log(covar_t$lagpopden+1)
  x%v%"laggdpgrowth"<-covar_t$lagGDPgrowth
  x
})

test<-ergm(tmp[[1]]~triadcensus(0:3)+gwdegree(.5,fixed=TRUE),estimate='MPLE')
test2<-ergm(tmp[[1]]~triadcensus(0:3)+degree(2:10),estimate='MPLE')
test3<-ergmMPLE(tmp[[1]]~triadcensus(0:3)+gwdegree(.5,fixed=TRUE))
test4<-ergmMPLE(tmp[[1]]~triadcensus(3)+nodecov("g77china")+sociality(nodes=c(185,38))+gwdegree(.5,fixed=TRUE))
head(test4$predictor)

input_network <- tmp

## Top 10 Carbon admission countries
## https://www.ucsusa.org/resources/each-countrys-share-co2-emissions
top10<-which(covar$country%in%c("China","United States of America","Russia","Japan",
                                "Germany","South Korea", "Iran", "Canada", "Saudi Arabia"))
data.frame(country=c("China","United States of America","Russia","Japan",
                     "Germany","South Korea", "Iran", "Canada", "Saudi Arabia"),top10)

top5<-which(covar$country%in%c("China","United States of America","Russia","Japan",
                               "Germany"))
top4<-which(covar$country%in%c("China","United States of America","Russia","Japan"))

top3<-which(covar$country%in%c("China","United States of America","Russia"))

top2<-which(covar$country%in%c("China","United States of America"))

top1<-which(covar$country%in%c("China"))


data.frame(country=c("China","United States of America","Russia","Japan",
                     "Germany"),top5)

# nodecov("laggdpgrowth")
#nodecov("lagpopden")
#nodecov("g77china")
#edgecov("log_d")
#gwdegree(.5,fixed=TRUE))
test4<-ergmMPLE(tmp[[1]]~triadcensus(3)+sociality(nodes=top5)+
                  nodematch("g77china")+
                  isolates)

head(test4$predictor)

model.formula <- net~triadcensus(3)+sociality(nodes=top5)+
  nodematch("g77china")+isolates;
model.terms <- colnames(test4$predictor)

input_network<-tmp
graph_mode <- 'graph'
group <- 'eu'
alpha.glmnet <- 1
directed <- TRUE
lambda <- NA
exvar <- NA

maxlag = 1
cdim <- length(model.terms)
lagmat = matrix(0,nc=cdim,nr=maxlag+1)
lagmat[1,]<-0
lagmat[2,]<-1
lagmat[1,grep("edgecov",model.terms)]<-1
lagmat[2,grep("edgecov",model.terms)]<-0
lagmat[1,grep("nodecov",model.terms)]<-1
lagmat[2,grep("nodecov",model.terms)]<-0


ylag = rep(1,maxlag)
out <- suppressWarnings(paramEdge(input_network,
                                  model.terms,
                                  model.formula,
                                  graph_mode='graph',
                                  group=group,
                                  intercept = NA, 
                                  exvar=NA,
                                  maxlag = maxlag,
                                  lagmat = lagmat,
                                  ylag = ylag,
                                  lambda = lambda, 
                                  method='bayesglm',
                                  alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+sociality(nodes=top4)+
  nodematch("g77china")+isolates;
model.terms2<-model.terms[1:5]
model.terms2<-append(model.terms2, model.terms[7:8])
lagmat2<-lagmat[,1:7]


out2 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+sociality(nodes=top3)+
  nodematch("g77china")+isolates;
model.terms2<-model.terms[1:4]
model.terms2<-append(model.terms2, model.terms[7:8])
lagmat2<-lagmat[,1:6]

out3 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+sociality(nodes=top2)+
  nodematch("g77china")+isolates;
model.terms2<-model.terms[1:3]
model.terms2<-append(model.terms2, model.terms[7:8])
lagmat2<-lagmat[,1:5]

out4 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+sociality(nodes=top1)+
  nodematch("g77china")+isolates;
model.terms2<-model.terms[1:2]
model.terms2<-append(model.terms2, model.terms[7:8])
lagmat2<-lagmat[,1:4]

out5 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+sociality(nodes=top5)+isolates;
model.terms2<-model.terms[1:6]
model.terms2<-append(model.terms2, model.terms[8])
lagmat2<-lagmat[,1:7]

out6 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+
  nodematch("g77china")+isolates;
model.terms2<-model.terms[1]
model.terms2<-append(model.terms2, model.terms[7:8])
lagmat2<-lagmat[,1:3]

out7 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

model.formula2 <- net~triadcensus(3)+sociality(nodes=top5)+
  nodematch("g77china");
model.terms2<-model.terms[1:6]
model.terms2<-append(model.terms2, model.terms[7])
lagmat2<-lagmat[,1:7]

out8 <- suppressWarnings(paramEdge(input_network,
                                   model.terms2,
                                   model.formula2,
                                   graph_mode='graph',
                                   group=group,
                                   intercept = NA, 
                                   exvar=NA,
                                   maxlag = maxlag,
                                   lagmat = lagmat2,
                                   ylag = ylag,
                                   lambda = lambda, 
                                   method='bayesglm',
                                   alpha.glmnet=1))

#Prints the coefficients
out$coef$fit
out2$coef$fit
out3$coef$fit
out4$coef$fit
out5$coef$fit
out6$coef$fit
out7$coef$fit
out8$coef$fit

#Creates a Coefficient Table. NOTE Model #'s don't match published version********
coef_list <- lapply(list(out, out2, out3, out4, out5, out6, out7, out8), function(out) {
  coefs <- coef(out$coef$fit) # Extract coefficients
  return(coefs)
})

# Get all unique coefficient names across all models
all_coefficients <- unique(unlist(lapply(coef_list, names)))

# Align all coefficients into a single dataframe
aligned_df <- sapply(coef_list, function(model_coefs) {
  # Match coefficients to the full list of all coefficients, filling missing ones with NA
  aligned <- setNames(rep(NA, length(all_coefficients)), all_coefficients)
  aligned[names(model_coefs)] <- model_coefs
  return(aligned)
})

# Convert to a proper dataframe
aligned_df <- as.data.frame(aligned_df)
colnames(aligned_df) <- paste0("Model", 1:8) # Name the columns
rownames(aligned_df) <- all_coefficients # Set rownames to coefficient names

write.csv(aligned_df, "Model_Comparison_Covars.csv")

#Then same thing for SE's
library(arm)
coef_list <- lapply(list(out, out2, out3, out4, out5, out6, out7, out8), function(out) {
  coefs <- se.coef(out$coef$fit) # Extract coefficient SE's
  return(coefs)
})

# Get all unique coefficient names across all models
all_coefficients <- unique(unlist(lapply(coef_list, names)))

# Align all coefficients into a single dataframe
aligned_df <- sapply(coef_list, function(model_coefs) {
  # Match coefficients to the full list of all coefficients, filling missing ones with NA
  aligned <- setNames(rep(NA, length(all_coefficients)), all_coefficients)
  aligned[names(model_coefs)] <- model_coefs
  return(aligned)
})

# Convert to a proper dataframe
aligned_df <- as.data.frame(aligned_df)
colnames(aligned_df) <- paste0("Model", 1:8) # Name the columns
rownames(aligned_df) <- all_coefficients # Set rownames to coefficient names

write.csv(aligned_df, "Model_Comparison_SEs.csv")
